

      real  function hs(phi,cutoff)
      real  phi,cutoff

      if (phi.ge.cutoff) then
        hs=1.0
      elseif (phi.le.-cutoff) then
        hs=0.0
      else
        hs=phi/( 2.0 *cutoff)+sin( 3.14159265358979 *phi/cutoff)/( 
     $    2.0 * 3.14159265358979 )+ 0.5 
      endif
      return
      end

      real  function rhoh(phi,cutoff)
      real  phi,cutoff,temp_hs,hs

      temp_hs = hs(phi,cutoff)
      rhoh=1000.0*temp_hs+1.225*(    1.0 -temp_hs)

      return
      end
